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Abstract 


In open canopies, the within-canopy flux from the ground surface and 
understory can account for a significant fraction of the total flux above the 
canopy. This study incorporates the important influence of within-canopy 
stability on turbulent mixing and subcanopy fluxes into a first-order clo- 
sure scheme. Toward this goal, we analyze within-canopy eddy-correlation 
data from the old aspen site in the Boreal Ecosystem - Atmosphere Study 
(BOREAS) and a mature ponderosa pine site in Central Oregon, USA. A 
formulation of within-canopy transport is framed in terms of a stability- 
dependent mixing length, which approaches Monin-Obukhov similarity the- 
ory above the canopy roughness sublayer. 

The new simple formulation is an improvement upon the usual neglect of 
the influence of within-canopy stability in simple models. However, frequent 
well-defined cold air drainage within the pine subcanopy inversion reduces 
the utility of simple models for nocturnal transport. Other shortcomings of 
the formulation are discussed. 


Key words: Canopy mixing, Cold air drainage, Similarity theory, Sub- 
canopy turbulence 

1 Introduction 

The need to understand and quantify biosphere-atmosphere exchange has led 
to substantial interest in within-canopy turbulence. Coherent eddies encom- 
passing the entire canopy often dominate the vertical transport and can lead 
to fluxes counter to the local vertical gradient within the canopy (Raupach 
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et al., 1996, Blanken et al., 1998, Finnigan, 2000). The primary source of 
canopy turbulence is often thought to be downward transport of turbulence 
from above the canopy. 

Modelling within canopy transport normally assumes a Lagrangian or 
Eulerian approach. This study follows the Eulerian approach. A common 
modelling approach extrapolates the mixing coefficient from the surface layer 
above the canopy downward through the canopy using a specified functional 
dependence on height (e.g., Shuttleworth and Wallace, 1985; Bonan, 1996). 
Modelling the flux at a given level includes first-order closure (e.g., Wilson et 
al., 1998), K-epsilon (Katul et al., 2004), second-order closure (e.g., Poggi et 
al., 2004) and large-eddy models (e.g., Shaw and Shumann, 1992, Albertson 
et al., 2001). These models generally require a flux-gradient approximation 
at some level of parameterization, often expressed in terms of a mixing length 
(e.g., Raupach et al., 1996; Wilson et al., 1998; Katul et al., 2004; Poggi et 
al., 2004). 

Using a mixing-layer analogy, Raupach et al. (1996) introduced the shear- 
length for describing turbulence above the canopy top. This length scale has 
been incorporated into several mixing-length formulations. For example, Wil- 
son et al. (1998) proposed a height-dependent mixing length that approaches 
the shear length scale near canopy top. Poggi et al. (2004) employed a length 
scale related to element diameter and the distance from the ground in addi- 
tion to the shear-length scale in the mixing length formulation. On the other 
hand, Massman and Weil (1999) focused on the influence of local leaf area 
index on the mixing length and related the mixing length inversely to the 
local drag coefficient such that the mixing length was smaller at levels with 
more local leaf area index. 

Less attention has been given to the mixing length for momentum and 
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heat within open canopies and the potential role of surface heating and cool- 
ing at the ground/understory surface. Buoyancy effects in the subcanopy 
have been included in a few LES models (e.g. Dwyer et al., 1997; Albertson 
et al., 2001). Brunet and Irvine (2000) investigated the influence of sta- 
bility on the shear-length scale for an open canopy and suggested that the 
primary effect of atmospheric stability is to modify the shear-length scale 
through change of wind speed and wind shear at the canopy top. However, 
they did not examine the influence of local stability on mixing within the 
canopy. Based on eddy-correlation data, Mahrt et al. (2000) found that the 
drag coefficient in the subcanopy at the old aspen site in BOREAS decreased 
with increasing subcanopy stability but did not increase from near-neutral 
to unstable conditions. However, they did not provide a formulation for such 
transport. 

Much of the work above is influenced by observations of the important 
contribution of coherent structures to the within-canopy flux, appearing to 
originate above the canopy, often conceptualized in terms of within-canopy 
sweeps and ejection (e.g. Paw U et al. 1992). Such structures do not pen- 
etrate the strongly stratified subcanopy associated with a open overstrory 
and clear nocturnal conditions (Mahrt, et al . 2000). In open canopies, day- 
time subcanopy turbulence may be generated more by local buoyancy than 
downward transport of turbulence from above the canopy. Resulting under- 
estimation of subcanopy turbulence can significantly influence estimates of 
land-surface exchange in open canopies. 

In this study, we analyze eddy-correlation data from two open canopies 
and construct a simple subcanopy formulation of the mixing length, which 
includes the influence of subcanopy stability. 
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2 Data 


This study analyses eddy-correlation data collected from a mature ponderosa 
pine site in central Oregon, USA (Schwartz et al., 2004). Turbulent fluxes 
were measured above the canopy and at two levels within the canopy, one 
in the crown space (10 m) and one in the trunk space (3 m). The flux at 
3m was measured at a subcanopy tower separated from main tower by 40 
m to avoid disturbance of the flow at the tower base. Wind speed and air 
temperature were measured at 5 levels (3 m, 6 m, 10 m, 20 m, 30 m) on the 
main tower. Air temperature was also measured on the subcanopy tower at 
1 m, 2 m and 3 m. 

This study also analyzes tower data collected at five levels on the main 
tower in the old aspen site in BOREAS (Sellers et al., 1995; Blanken et 
al., 1997). Momentum fluxes were computed from 30-min records. The 
instrumentation at the old aspen site is detailed in Blanken et al. (1997). 
Here we include only the post leafout period when the data were the most 
complete. 

Both canopies are characterized by a LAI of about 3 (Table 1), which we 
refer to as open canopies. While such values are not low compared to sparser 
canopies with a LAI of unity or less. Here, “open” refers to large diurnal 
variation of atmospheric stability within the canopy layer due to radiative 
heating and cooling of the ground/understory surface. The nocturnal tem- 
perature profiles at the old aspen site correspond to a strong surface inversion 
in the lowest 5 m and weaker stratification above (Mahrt et al., 2000). The 
stratification of the daytime aspen canopy layer is near neutral. The diurnal 
variation of the stratification at the pine site is much greater than that at 
the aspen site in spite of similar LAI values for the overstory. The strong 
diurnal variation of stability at the pine site is due partly to more clumping, 
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dryer surface conditions, less understory, higher sun angles, much less cloud 
cover and minimal integrable atmospheric water vapour content, all acting to 
increase the radiative heating and cooling of the ground surface/ understory. 

Fluxes are computed using deviations from 10-minute averages and then 
are averaged over one hour to reduce random flux errors. Fluxes in the sub- 
canopy are subject to some uncertainty from possible flux loss due to path- 
length averaging and expected horizontal heterogeneity above the nonuniform 
understory. Based on our assessment of momentum flux loss due to path- 
length averaging, by comparing with three-dimensional hotfilm anemometry 
(unpublished), we have determined that pathlength average can lead to sig- 
nificant momentum flux loss in the lowest meter above short grass and that 
existing correction formulas are not suitable. However, we have not carried 
out suitable measurements in the subcanopy. 


Table 1 Site description. 


Sites 

Average canopy 
height (m) 

LAI 

displacement 
height (m) 

roughness 
length (m) 

Old aspen 

20.1 

3.0 

13.4 

2.0 

Pine 

15.5 

3.3 

11 

1.2 


The displacement height for the aspen site (Table 1) is taken from Naka- 
mura and Mahrt (2001) while the displacement height for the pine site was 
estimated by trial and error to produce the best fit to Monin-Obukhov sim- 
ilarity theory at the top of the tower for near-neutral conditions. 


6 


53 



3 Transport in the canopy layer 

3.1 Mixing length 

We express the within-canopy flux in terms of an eddy diffusivity 

W = -K m ^, ( 1 ) 

In the data analysis below, the momentum flux will be computed as the 
magnitude of the vector momentum flux and the mean shear will be computed 
as the magnitude of the vector shear without concern for alignment of the 
wind and shear vectors. 

The eddy diffusivity is formulated as 


— lm.U* 


( 2 ) 


where l m is the mixing length for momentum and n* is the height dependent 
friction velocity in the canopy defined as the square root of the magnitude 
of the momentum flux. The mixing length for momentum 

u, 


< 3 > 

is evaluated directly from eddy-correlation and wind profile measurements at 
the pine and aspen sites. 

The mixing length for momentum at the aspen site was computed from 
the observed fluxes and vertical gradients for three different stability classes, 
based on the stability at the top of the tower (39- m level): the unstable case 
(—1 < (z — d)/L < —0.1), the near neutral case (—0.1 < (z — d)/L < 0.1) 
and the stable case (0.1 < (z — d)/L < 1). For this study, transition periods 
are eliminated and we analyze only daytime data between 1000 and 1600 
local time and nocturnal data between 2100 and 0500 local time. Vertical 
gradients are computed from simple finite differencing. 
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The resulting mixing length for momentum, averaged over all of the quali- 
fying records, increases slowly with height in the subcanopy and more rapidly 
with height above the canopy (Figure 1). The mixing length is much smaller 
in the subcanopy in stable conditions compared to neutral and unstable con- 
ditions, probably due to weaker downward transport of turbulence energy 
into the canopy layer. The profile of the mixing length for the pine site is 
based on only three levels and is not shown. 

To formulate the height dependence of the mixing length, we recognize 
that height above the ground surface influences the within-canopy mixing 
while the mixing length must approach Monin-Obukhov similarity theory in 
the surface layer above the canopy roughness sublayer. For neutral conditions 
within the canopy, we propose the simple height dependence 

^mn " (3z (4) 

where l mn is the neutral value of the mixing length for momentum, z is 

the height above the ground surface. (3 is an empirical coefficient to-be- 

determined and not necessarily equal to the von Karman constant. The 
nondimensional coefficient (3 is estimated from the observed profiles, yielding 
0.18 at the aspen site and 0.30 at the pine site. Either f3 depends on canopy 
architecture, or, it cannot be estimated from our observations within a factor 
of two. 

In order that the mixing length satisfy Monin-Obukhov similarity theory 
in surface layer above the canopy, (3 at the top of the roughness sublayer, z r , 
should approach 

(3 = 0A(z r - d)/z r (5) 

where d is the displacement height and z r can be estimated following Raupach 


8 


55 



• ■ ■ unstable 
neutral 
■ ■ stable 


Figure 1: Height dependence of the mixing length for momentum estimated 
from vertical profiles at the old aspen site for the post leafout period. 
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(1994) as z r = d+2(h — d) where h is the canopy height. Assigning d = C h, 

(3 = 0.4[1 - t^q\ (6) 

As an exercise, we specify the displacement height to be 2/3h (C=2/3), in 
which case (3 becomes 0.2, which is closer to value at the aspen site. For 
the open canopies of this study, the value of C may be smaller leading to a 
somewhat larger value of (3. Although there is no reason to expect the profile 
fit to Eq. 4 to agree with Eq. 6, we note that they both predict comparable 
magnitudes. Below, we use values of /3 from the profile fit since it avoids 
estimating the height of the roughness sublayer and displacement height. 

3.2 Stability influence on momentum transfer 

In analogy to Monin-Obukhov similarity theory, the stability-dependent mix- 
ing length can be expressed as 


7 - 

'"171 ~~ , ■ 

0m 

Here, the nondimensional wind shear is defined as 

Pzdu/dz 

4>m = • 

u* 

We do not expect (f> m to closely follow Monin-Obukhov similarity theory 
within the canopy because the stress varies rapidly with height. The flux- 
gradient relationship in the subcanopy will be limited in different ways by the 
distance from the ground surface/understory, the spacing between the trunks 
and the distance from the overstory. However, both the heat flux and the 
friction velocity are important influences so that the local Obukhov length, 
A, is one of the governing length scales. So that the stability approaches 
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Monin-Obukhov similarity theory at the top of the roughness sublayer, we 
formulate the subcanopy stability parameter as 


Z.Zr-d 

A 1 * r ' 


( 9 ) 


where the ratio in parentheses allows continuity of the local canopy stability 
functions to that above the roughness sublayer. 

The nondimensional wind shear (Figures 2-3) varies more slowly with sta- 
bility than predicted by the normal stability functions for Monin-Obukhov 
similarity theory at both the aspen and pine sites. The relationship between 
the nondimensional shear and stability exhibits substantial scatter, which 
could be due partly to the difficulty of measuring fluxes in the subcanopy 
(Section 2) and the omission of additional physical influences. In spite of 
the large scatter, inclusion of the influence of subcanopy stability is needed, 
although more data is required to determine the optimum stability function. 
The disagreement between the observations and the usual Monin-Obukhov 
stability functions is greatest at the lowest level of the aspen site, just above 
the understory, where the stability dependence is almost undetectable. Inclu- 
sion of a displacement height for the thick understory of 2-m height increases 
the rate at which the nondimensional shear changes with stability and thus 
increases the stability influence. However, this increase is not nearly large 
enough to allow reasonable approximation with the usual Monin-Obukhov 
stability functions. However, the agreement with the Monin-Obukhov sta- 
bility functions improves with height (Figure 2). 

The smaller observed nondimensional shear (greater mixing) for stable 
conditions, compared to the prediction by the usual stability functions, is 
presumably due to mechanical generation of mixing by the vegetation ele- 
ments, not included in Monin-Obukhov similarity. However, the nondimen- 
sional shear above the canopy is also observed to be smaller than predicted 
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by linear dependence on stability for strong stability (e.g. Beljaars and Holt- 
slag, 1991; Cheng and Brutsaert, 2005). The larger nondimensional shear 
for unstable conditions, most evident at the pine site, could be due to the 
partial constraint of the convective eddies by the height of the overstory, 
which enters as an additional length scale not included in Monin-Obukhov 
similarity. 


3.3 Stability influence on heat transfer 


The pine canopy is characterized by unstable stratification (Figure 4) and 
significant upward heat flux in the daytime and strong stable stratification 
at night. The diurnal variation of the within-canopy stability is less at the 
aspen site as discussed in Section 2. 

The heat transport is expressed in terms of the eddy diffusivity for heat 

ew = -K h ^ (io) 


The eddy diffusivity for heat can be related to the mixing length for heat. 
The mixing length for heat for near neutral conditions is difficult to estimate 
from data because of large scatter. We equate the neutral values of the 
mixing lengths for heat and momentum but allow them to have different 
stability functions. 

The nondimensional temperature gradient is defined as 


f3zd0/dz 


where 


ea- 


rn'd' 

u„ 


( 11 ) 


( 12 ) 


was evaluated at the two eddy-correlation levels at mature pine site, one 
in the crown space (10m) and the other in the trunk space (3m). This 
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Figure 2: The nondimensional shear as a function of stability for 3 different 
levels at the Aspen site and the nondimensional temperature gradient at 
the 18-m level. The top of the canopy is approximately 20 m. The solid 
line corresponds to the stability function from Dyer (1974) applied to 7 z/L, 
where z is the height above the ground surface. 7 = 


13 






| (a) 10m (Crown space) (b) 3m (Trunk space) 



Figure 3: The nondimensional shear and nondimensional temperature gra- 
dient for two levels at the pine site. The top of the canopy is approximately 
15 m. The solid line corresponds to the stability function from Dyer (1974) 
applied to 7 z/L. 7 = 
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relationship was also evaluated at the 18-m level at the aspen site where the 
heat flux and vertical temperature gradients are best defined for evaluation 
of the relationship. Evaluation of vertical gradients must recognize that the 
time-averaged temperature can vary horizontally in the subcanopy on the 
microscale, especially close to the ground surface at the pine site. We use 
only data where the horizontal temperature difference at 3m between two 
subcanopy towers and the main tower is less than 0.5 K, which eliminates 
some of the most stable cases. 

At the pine site, the nondimensional temperature gradient more closely 
follows the Monin-Obukhov stability function compared to that for momen- 
tum (Figure 3). Apparently, the influence of diabatic heating in the canopy 
layer and associated height-dependence of the heat flux does not strongly 
modify the dependence of the flux-gradient relationship on stability. How- 
ever, the influence of stability on the nondimensional temperature gradient 
is not as well defined for the aspen site (Figure 2). At both sites, the scatter 
is large for stable conditions, indicating large errors in the flux calculations 
and/or additional influences on the nondimensional temperature gradient, 
not included in z/A. 

Although the nondimensional temperature gradient for heat exhibits a 
stronger dependence on within-canopy stability compared to the nondimen- 
sional shear, the present data do not allow confident assessment of the ver- 
tical structure of the differences between the mixing lengths for heat and 
momentum. 

With the above findings and uncertainties, we formulate the eddy diffu- 
sivity for heat as 

K h = (13) 

0/1 

where l mn is the mixing length for momentum for neutral conditions (Eq. ??), 
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m* is the height-dependent friction velocity and <ph is the nondimensional tem- 
perature gradient with the usual Monin-Obuhkov stability functions, which 
formally enter here as the stability dependence of the mixing length. Since 
the dependence of the nondimensional shear on stability is less certain than 
that for temperature, we choose to follow the traditional top down approach 
of extrapolating the stress downward from the surface layer using a profile 
model. Ultimately, a separate equation for the momentum budget is needed 
for the canopy layer. 

4 Influence of LAI on canopy momentum flux 

Estimation of the diffusivity for heat from the stability-dependent mixing 
length in the previous section requires information on the mean profile of the 
local friction velocity in the canopy, which we now parameterize based on 
the observed behavior. The ratio of the subcanopy friction velocity to that 
above the canopy increases with increasing instability (Figure 5), presumably 
due to greater mixing in the daytime as suggested by larger within-canopy 
mixing lengths in the daytime (Figure 1). We will include the influence of 
the vertical structure of the LAI on the height dependence of the friction 
velocity. For the aspen site, the downward accumulated Leaf Area Index 
( LAI acc ) is estimated based on the canopy structure of the aspen (Blanken 
et al. 1997). We then relate the momentum flux in the mean wind direction 
to stability and canopy architecture as 

ul = ul 0 exp[-2a{LAI acc f h ] (14) 

where a is a stability-dependent coefficient evaluated below and u* 0 is the 
friction velocity in the surface layer above the canopy. The exponential form 
allows the momentum flux to approach the surface layer value above the 
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canopy. The square root dependence on the accumulated LAI was determined 
as a reasonable fit to near-neutral conditions. Although Eq. 14 adequately 
simulates the vertical profile of local friction velocity at the aspen site, the 
significantly deviations of the observed from the predicted values at the pine 
site (Table 2) either indicate inadequate vertical resolution of the flux profiles 
at the pine site or imply that the local friction velocity for near-neutral 
conditions are affected by factors other than accumulated LAI. The pine site 
is characterized by greater foliar clumping and clumping of trees, allowing 
easier vertical exchange and greater within-canopy momentum flux compared 
to Eq. 14. 


Table 2. Vertical profile of the height-dependent friction velocity nor- 
malized to the value above the canopy for neutral conditions. The ratio is 
defined as the local u * divided by friction velocity above the canopy, u* 0 . 


Aspen site 

Pine site 

z (m) 

LAI acc 

Ratio 

z(m) 

LAI acc 

Ratio 

observed 

estimated 

observed 

estimated 

28 

0 

i 

i 

30 

0 

i 

i 

18.6 

1. 

0.49 

0.37 

10 

1.5 

0.62 

0.29 

5.8 

3. 

0.19 

0.18 

3 

3 

0.13 

0.18 

0.45 

6.5 

0.05 

0.08 






Using the observed fluxes, we have computed a from Eq. 14. We have 
visually fit the stability dependence of a for unstable conditions as (Figure 

5 ) 


a (1-42/L) 0 - 25 ' ^ 

Application of this stability dependence in Eq. 14 leads to reasonable pre- 
diction of the ratio of the subcanopy stress to the stress above the canopy for 
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unstable conditions at both the aspen and pine sites. For stable conditions, 
a appears to be independent of stability (Figure 5). While the above inclu- 
sion of differences between stable and unstable conditions is rather crude, it 
is expected to be an improvement upon the usual practice of neglecting the 
influence of stability on the subcanopy transport. However, with the frequent 
occurrence of drainage flow and capping inversion at the pine site (Section 
5), the vertical communication is too weak to apply Eq. 14 and a separate 
time-dependent momentum equation is needed near the surface. 

5 Nocturnal structure of the pine subcanopy 

Part of the large scatter for the nocturnal flux-gradient relationship at the 
pine site is due to cold air drainage and wind directional shear. While some 
cold air drainage also occurred at the aspen site (Lee, 1998), it seemed to 
be less common and have less impact on the vertical structure of the flow 
within the canopy compared to the pine site. With cold air drainage at the 
pine site (westerly flow near the surface), the nocturnal stratification is much 
stronger and concentrated at the top of the cold air drainage, between 10 and 
20 m (Figure 4). For individual hours, the stratification is often confined to 
a single 5-m observational layer but is vertically spread in the composited 
profile of Figure 4. As a result, the nocturnal vertical temperature gradients 
at the 10 m level are probably often seriously underestimated. A network 
of wind and temperature measurements indicates cold air drainage down the 
eastward descending slope, just west of the tower site. In the absence of 
cold air drainage, the subcanopy flow is generally from the south and the 
stratification above the under story is weaker. 

The joint frequency distribution between the wind speed at 30 m and the 
temperature difference in canopy between 3m and 20m (Figure 6) indicates 
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39m 


(a) Pine 10 m (b) Pine 3 m 






Figure 5: The ratio of friction velocities as a function of stability where L a j c 
is the Obukhov length based on fluxes at the 39-m level, assumed to be in 
the surface layer. The solid line represents the prediction of the ratio based 
on the fit to Eq. 15 for unstable conditions and independent of the stability 
for stable conditions. 
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Figure 6: Joint frequency distribution of the 30 m wind speed and the canopy 
inversion strength between 3 and 20 m (DT) at the pine site. Color pattern 
is only for visualization. 
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Figure 7: Joint frequency distribution of subcanopy 3-m wind direction and 
the canopy inversion strength between 3 and 20 m (DT) at the pine site. 
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that most of the nights with inversions greater than a few degrees occur with 
weak winds above the canopy and that strong winds generally correspond to 
weaker inversions of less than a few degrees. Almost all of the strong inver- 
sions occur with westerly subcanopy drainage flow. Weak canopy inversions 
correspond to southerly subcanopy flow more aligned with the flow above the 
canopy (Figure 7). The systematic development of subcanopy drainage flows 
and associated directional shear and strong inversions at the pine site prob- 
ably account for much of the large scatter in the flux-gradient relationship. 
Inclusion of the directional shear associated with the drainage flow requires 
a separate momentum equation for the subcanopy. 

6 Summary 

The structure of turbulent transport within the canopy has been examined 
with eddy-correlation measurements from a relatively open ponderosa pine 
site and from an aspen site. The dependence of the nondimensional temper- 
ature gradient and mixing length for heat on stability approximately follow 
MO similarity theory even though the assumptions required for such simi- 
larity theory are not met. The nondimensional shear and mixing length for 
momentum deviate more from MO similarity. Deviations are most significant 
close to the top of the understory where the mixing coefficient depends much 
less on stability than would be estimated from the MO stability functions. 
Inclusion of the influence of subcanopy stability on subcanopy transfer is 
expected to improve formulation of transport of heat, momentum and pre- 
sumably trace gases such as carbon dioxide. 

Subcanopy drainage flows with strong capping inversions are common at 
the pine site and contribute to deviations from simple similarity theory. Ac- 
curate formulation of the within-canopy transport for these situations would 
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require a separate subcanopy momentum equation. We are collecting data 
with better vertical resolution to determine the influence of this inversion on 
decoupling between surface fluxes and fluxes above the canopy. 

In spite of these difficulties and uncertainties with the above analyses, 
it is possible to improve upon existing simple formulations of within-canopy 
mixing, which completely neglect the dependence of within-canopy mixing on 
the stability. For example, an eddy diffusivity for heat has been constructed 
(Eq. 13) from the mixing length for neutral conditions (Eqs. 4, 5), the 
usual Monin-Obukhov stability functions and the formulation of the profile 
of the within-canopy friction velocity scale (Eqs. 14, 15), which depends on 
downward accumulated LAI and within-canopy stability. 

More complete subcanopy flux profiles and their microscale variability 
are needed to provide more dependable profiles of the mixing length for heat 
and the dependence of the nondimensional shear on canopy architecture and 
within-canopy stability. 
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